Quantum chaos and thermalization in gapped systems 
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We investigate the onset of thermalization and quantum chaos in finite one-dimensional gapped systems 
of hard-core bosons. Integrability in these systems is broken by next-nearest-neighbor repulsive interactions, 
which also generate a superfluid to insulator transition. By employing full exact diagonalization, we study chaos 
indicators and few-body observables. We show that with increasing system size, chaotic behavior is seen over 
a broader range of parameters and, in particular, deeper into the insulating phase. Concomitantly, we observe 
that, as the system size increases, the eigenstate thermalization hypothesis extends its range of validity inside 
the insulating phase and is accompanied by the thermalization of the system. 
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The use of ultracold quantum gases to realize strongly cor- 
related phases of matter has been at the forefront of research 
during the last decade J2]. Such systems will not only help 
to understand phases and models introduced in the past, but 
also will create and will provide the means to investigate new 
exotic phases. Another major promise of the field is the pos- 
sibility for studying the dynamics of correlated quantum sys- 
tems in very controlled ways. For short times, the collapse 
and revival of phase coherence was analyzed in Ref. yj] af- 
ter an interaction quench from the Bose-Einstein condensate 
to the Mott insulating regime. Ab initio calculations of var- 
ious bosonic and fermionic lattice models reproduced those 
observations iH-Ot]. A natural question that follows, given the 
fact that these systems are nearly isolated, is whether conven- 
tional statistical ensembles can describe experimental observ- 
ables after relaxation. 

Recent experiments have found that relaxation toward ther- 
mal equilibrium takes place in certain setups, but not in others 
@]. The lack of thermalization in the latter case can be asso- 
ciated with the proximity to integrability [7]. At integrability, 
several works have shown that thermalization is not expected 
to occur JH], except for special conditions and/or observables 
JUI^]. For generic nonintegrable systems, thermalization is 
expected to occur and follows from the eigenstate thermaliza- 
tion hypothesis (ETH) JToLITTIl . 

Interestingly, in numerical calculations for quenches across 
a gapless to gapped phase, thermalization did not occur in the 
gapped side, even when the system was nonintegrable 0]. As 
a matter of fact, ETH, which was shown to hold in several 
nonintegrable gapless systems 10, [TUl . has been questioned 
for gapped systems lfl2l [l3ll . Questions raised include the 
proximity to the atomic limit, finite-size effects Il2tl . and the 
effects of rare states II 1 3(1 within the insulating phase. Other 
studies that deal with spin and fermionic systems have found 
that relaxation toward equilibrium occurs faster close to a crit- 
ical point J3. The fact that, away from the ground state, 
these systems are, in general, not insulating, even if gaps are 
present in the spectrum, renders the debate more interesting 
still. Why then would generic nonintegrable gapped systems 
behave differently from gapless ones? Remarkably, in disor- 
dered systems, insulating behavior may take place away from 
the ground state [ Oil, and those ones could indeed behave dif- 
ferently. 



In this Rapid Communication, we use various measures, 
such as quantum chaos indicators lfl6ll and eigenstate expec- 
tation values of experimental observables 131 to understand 
whether thermalization should occur in gapped systems. This 
question is of interest for current experiments with ultracold 
gases, where systems with insulating ground states such as the 
bosonic and fermionic Mott insulators are studied. We show 
that thermalization does occur in the gapped side of the phase 
diagram and that, as the system size increases, thermalization 
is observed deeper into that phase. We also find that ETH 
holds for systems with larger gaps, as the system size is in- 
creased. This supports the view that thermalization occurs in 
generic nonintegrable systems independent of the presence or 
absence of gaps in the spectrum. One does need to be care- 
ful with temperature effects and finite-size effects, which may 
be more relevant close to integrable points 13] and in systems 
with gaps lfl2tl . We also study the long-time dynamics of these 
gapped systems and address its universality. 

We focus our study on the one-dimensional hard-core 
boson (HCB) model with nearest-neighbor hopping t, and 
nearest- and next-nearest-neighbor (NNN) interaction V and 
V , respectively. The Hamiltonian is written as 

L 

flt = £){-*(S&+i+H.c.) (1) 

i=l 

where standard notation has been used 13]. We restrict our 
analysis to lattices with 1/3 filling (Nb — £/3); t — 1 sets 
the energy scale, V = 6, and V' is varied (0 < V' < 9). 
For V' < VJ. = 3, the ground state is a gapless superfluid, 
whereas for V' > V' c = 3, it is a gapped insulator IU7I1 . 

Eigenstate expectation values (EEVs) of different observ- 
ables, as well as the nonequilibrium dynamics and thermody- 
namics of these systems, are determined using full exact di- 
agonalization of the Hamiltonian in Eq. (Q~|i. We study lattices 
with up to 24 sites and 8 HCBs, which correspond to a total 
Hilbert space of dimension D = 735 471. We take advantage 
of the translational symmetry of the lattice to independently 
diagonalize each Hamiltonian block with total momentum k; 
the largest block in k space has dimension Dk = 30 667 lfl8ll . 
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The Hamiltonian in Eq. (Q~|i is integrable when V' = 0, 
whereas the addition of NNN interaction leads to the onset of 
chaos. Note that the integrable-chaos transition may occur al- 
though random elements are nonexistent in the Hamiltonian 
lfl6ll . We start our study by addressing how quantum chaos in- 
dicators, such as level spacing distribution, level number vari- 
ance, and inverse participation ratio (IPR) HH, change as one 
moves away from the integrable point and eventually enter the 
gapped region by increasing V '. The outcomes are shown to 
support our results on thermalization. 

Level spacing distribution and level number variance are 
obtained from the unfolded spectrum of each k sector sepa- 
rately. The first is a measure of short-range correlations, and 
the second is a measure of long-range correlations lfl9ll . For 
integrable systems, the distribution of spacings s of neigh- 
boring energy levels may cross, and the distribution is Pois- 
sonian Pp(s) = cxp(— s), while for nonintegrable systems, 
level repulsion leads to the Wigner-Dyson distribution. The 
form of the latter depends on the symmetries of the system. 
Here, it coincides with that of ensembles of random matrices 
with time reversal invariance, the so-called Gaussian orthog- 
onal ensembles (GOEs): Pgoe(s) = (tts/2) exp(— 7rs 2 /4). 
The level number variance is defined as X 2 (7) = (N(l, e) 2 } — 
(N(l, e)) 2 , where N (I, e) is the number of states in the energy 
interval [e, e+l] and (.} is the average over different initial val- 
ues of e. For a Poisson distribution Sp(/) = /, while for GOEs 



where b 



13+2 
0+1 



in the limit of large I, £ GOE (0 = 2[ln(27rZ)+7+l-7r 2 /8]/7r 2 



where 7 is the Euler constant. 




FIG. 1 : (Color online) (a) Parameter /3 of the Brody distribution used 
to fit the average of the level spacing distributions over all sectors 
with k — 1, . . . , \_L/2\ . The inset shows the gap (between the lowest 
energy states) times L vs V'. (b) Level number variance averaged 
over the same k sectors, for L = 24, and compared to the Poisson 
and GOE results, (c) Inverse participation ratio in the mean-field 
(mf) (top) and momentum (bottom) basis vs energy per site; L = 24, 
k = 2, and D k = 30 624. 

Figures [Tfa) and [Tib) show results for P(s) and £ 2 (0> re- 
spectively, for different values of V'. The level spacing distri- 
bution is parametrized by /?, which is used to fit P(s) with the 



Brody distribution J2fj, P B (s) = (fi + l)&s /3 exp (-6s' 3+1 ), 



0+1 j— . 

II18H . Two transitions are veri- 
fied: (i) from integrable [j3 -> 0, S 2 (0 — > I] to chaotic 
[/3 -> 1,S 2 (0 Egoe(01 as v ' increases from V = 0, 
and (ii) a departure from chaoticity for large values of V'. 
Figure [JJ a) also shows a strong dependence of the results on 
the system size. For larger systems, smaller values of V suf- 
fice for the first integrable-chaos transition, and larger values 
of V' are required for the second transition. It is an open 
question whether, in the thermodynamic limit, any value of 
V 7^ would implicate chaoticity for these systems. How- 
ever, for our finite systems, we do find an overlap between the 
gapped phase and the chaotic regime. The behavior of the gap 
A (times L) vs V' is depicted in the inset in Fig. 01 a )- No- 
tice the kink around V' c ~ 3, which signals the onset of the 
superfluid-insulator transition. 

The IPR measures the level of derealization of the eigen- 
states |2lll . Contrary to the two previous quantities, IPR de- 
pends on the basis chosen for the analysis. For an eigenstate 
\ip a ) of Eq. ([TJ written in the basis vectors \<fij) as \ip a ) = 
Ef=i wehaveIPR Q (£f = \ ^J 4 )" 1 . InFig.[lJc), 

we investigate IPR in two bases: the mf basis (IPR m f), where 
the are the eigenstates of the integrable Hamiltonian 

(V' = 0); and the k basis (IPRfe), where the |0j)'s are the total 
momentum basis vectors. Small vs larg e values of IPR m f sep- 
arate regular from chaotic behavior j2211 and signal derealiza- 
tion during the first integrable-chaos transition. The reduction 
of IPRfe indicates localization in k space, which in our model, 
results from approaching the atomic limit in the presence of 
translational symmetry. Hence, it explains the departure from 
chaoticity for V > V' c HH. 

The structures of the eigenstates are intimately connected 
to the thermalization process. As stated by the ETH, ther- 
malization is expected to occur when the expectation values 
of experimental observables (with respect to eigenstates of 
the Hamiltonian that are close in energy) are very similar to 
each other and, hence, are equal to the microcanonical aver- 
age, that is, thermalization occurs at the level of eigenstates. 
This is certainly the case with GOEs, where the amplitudes c 3 a 
for all |^ Q )'s are independent random numbers. GOE eigen- 
states in any basis then lead to IPRgoe ~ -Dfe/3 l22ll . For our 
lattices, IPR m f and IPR^ may approach the GOE value only 
in the middle of the spectrum, as expected for systems with 
finite-range interactions lfl6ll23ll . States at the edges, even in 
the chaotic limit, are more localized. In addition, as shown in 
the panels of Fig.[Tfc), the values of IPRs for eigenstates close 
in energy fluctuate significantly as one moves away from the 
chaotic limit, namely, as V' — s- 0, where the system becomes 
localized in the mf basis, and for V ^> V' c , where the system 
becomes localized in the k basis. Thus, ETH is not expected 
to hold in these regions. On the contrary, for intermediate 
values of V', IPR m f and IPRfe become smooth functions of 
energy, especially for E < 0. We may, therefore, anticipate 
compliance with the ETH even after the opening of the gap. 

To verify the relation between ETH and the chaos measures 
analyzed previously, we have studied the EEVs of one- and 
two-body observables for different values of V' as one crosses 
the superfluid to insulator transition. We have found qualita- 
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tively similar results for them, so we only report here on the 
kinetic energy [K = ^\ —t (b\bi + i + VL.c^j] and the mo- 
mentum distribution function [n(k)]. Both are one-body ob- 
servables, although the former one is local, while the latter is 
not. K and n(k) are routinely measured in cold gases experi- 
ments. 




FIG. 2: (Color online) (a) EEVs of K (top) and n(k) (bottom) vs 
energy per site for the full spectrum (which include all momentum 
sectors). Results are shown for four different values of V' and L = 
24. Panels (b) and (c) average relative deviation of the EEVs of (b) 
K and (c) n(k = 0) with respect to the microcanonical result vs V' , 
for T = 3 (see text) and three different lattice sizes. T = 3 (L = 24) 
corresponds to E = -10.7 for V = 1, E = -10.0 for V' = 3, 
E = -11.7 for V = 5, and E = -19.3 for V = 9. 

Figure |2ja) depicts the EEVs of K and n{k) for all eigen- 
states of the Hamiltonian and for different values of V'. For 
small values of V' (V < 2 for L = 24), there are large 
fluctuations of the EEVs of both observables over the entire 
spectrum. As V increases and one departs from integrability, 
these fluctuations reduce in the center of the spectrum, and 
ETH becomes valid. To increase V even further (beyond the 
superfiuid to insulator transition), increases the fluctuations of 
the EEVs once again as the eigenstates begin to localize in 
k space. These results are in clear agreement with what we 
expected based on the chaos indicators lfl8tl . 

To be more quantitative, we study the average deviation of 
the EEVs with respect to the microcanonical result (A mlc ). 
For an observable O, we define A mlc O = (Y^ a \Oaa — 
Omic\)/(J2 a Oaa)> where the sum is performed over the mi- 
crocanonical window and O aa are the EEVs of O. The mi- 
crocanonical expectation values are computed as usual. We 
average over all eigenstates (from all momentum sectors) that 
lie within a window [E - AE, E + AE], and take AE = 0.1. 
We have checked that our results are independent of the ex- 
act value of AE in the neighborhood of AE = 0.1. Here, 
we select E such that, for different values of V', the effective 
temperature (T = 3) is the same for all systems sizes l24ll . 

Results for A mlc K and A mlc n{k = 0) are presented in 
Figs. |2|b) and [2c) lfl8ll . They show that, in general, as the 
system size increases: (i) The average deviations of both ob- 



servables decrease, and (ii) the upturn that occurs as localiza- 
tion starts to set in k space moves toward larger values of V' . 
A comparison between the lower and upper panels in Fig. [2] 
also shows that where A mlc K and A mlc n(k = 0) are mini- 
mal, so are the maximal fluctuations of K and n(k = 0) in the 
individual eig enstates, and they decrease with increasing sys- 
tems size 1181 . Hence, ETH is valid in that regime and we find 
no evidence of rare states lfl3ll . The above results are in agree- 
ment with the chaos measure predictions and indicate that, for 
thermodynamic systems, ETH may be valid, away from the 
edges of the spectrum, even if one is deep into the insulating 
side of the phase diagram. 

Equipped with this knowledge, we are now ready to study 
the dynamics of such systems after a quench. Our initial states 
are always selected from the eigenstates of ([T]i with t = 1, 
V = 6, V( ni and zero total momentum, and then we quench 
V( ni — > V'. After a systematic analysis, we have found that 
the short-time dynamics depends strongly on the initial state 
and the final Hamiltonian, so we will focus here on the long- 
time dynamics and the outcome of relaxation. As discussed 
previously JJ, 0, [Ol, after relaxation, observables are well 
described by the diagonal ensemble Odiag = J2 a \C a \ 2 O aa , 
where C a is the overlap of the initial state with eigenstate a 
of the Hamiltonian. 
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FIG. 3: (Color online) Dynamics of the normalized difference be- 
tween the evolving expectation values of K (left panels) and n(k) 
(right panels) and the diagonal ensemble prediction. Average over 
the evolution of nine initial states selected from the eigenstates of 
the Hamiltonian with V{ ni = 0, 1, . . . , 9 (excluding the V used for 
the dynamics). The nine states for each V' were chosen such that 
the (conserved) energies during the time evolution are the same in all 
cases, and T = 3. Given the energy of the initial state in the final 
Hamiltonian E — (ipini\H\ipi n i), T is computed by following Ref. 
& 

In Fig. [3] we show the normalized difference between 
the time-evolving expectation value of K and n(k) and 
the diagonal ensemble prediction SK (left panels) and 
6nk (right panels), respectively. We define 5K(t) = 
\K(t) - K diag \/\K diag \ and <5n fe (r) = (J2k\ n ( k ' T ) ~ 
n<Hag(k)\)/ (J2k n diag{k)). In order to verify the universal- 
ity of our results, for each value of V used for the dynamics, 
we prepared nine initial states selected from the eigenstates of 
the Hamiltonian with different values of V( ni (excluding V') 
and studied the dynamics for all of them. The long time dy- 
namics was found to be very similar, independent of the initial 
state. In Fig. [3] we depict the average over those nine different 
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time evolutions fl8ll . These plots show that after long times, 
observables relax to values similar to those predicted by the 
diagonal ensemble, and those predictions become more accu- 
rate with increasing system size. Only for V = 9 do we find 
large time fluctuations of K, which is a consequence of the 
approach to localization in k space. However, even these time 
fluctuations decrease with increasing system size. 

Once it is known that the relaxation dynamics brings ob- 
servables to the values predicted by the diagonal ensemble, 
and that the accuracy of that prediction improves with increas- 
ing system size, all we need to do to check whether the system 
thermalizes or not is to compare the predictions of the diag- 
onal ensemble with the microcanonical ones. For larger sys- 
tems, one could also compare with the canonical ensemble, 
but this is not adequate here due to finite-size effects 10]. 




FIG. 4: (Color online) Relative differences between the predictions 
of the microcanonical and diagonal ensembles for (a) K and (b) n(k) 
vs V', for T = 3.0 and T = 5.0. Results are shown for L = 21, 
Nb = 7 and L = 24, ATj = 8. The diagonal ensembles correspond 
to the quenches from the nine V( ni of Fig.[3]to V' , and the curves are 
averages over the nine relative differences. Relative differences are 
computed in exactly the same way as 8K and Sn^ in Fig. [3] 

In Fig. |U we depict the comparison between diagonal and 
microcanonical ensembles. Results are shown for the same 
set of quenches and initial states presented in Fig. [3] and for 
an additional set of initial states such that the effective tem- 
perature of the relaxed systems is a bit higher, namely, T = 5 



lfl8ll . The results for T — 3 and T = 5 are in qualitative agree- 
ment with each other for our two observables of interest. They 
show that the microcanonical ensemble predicts the outcome 
of the relaxation dynamics with high accuracy for the inter- 
mediate values of V, where ETH was shown to be valid (cf. 
Fig. |2j and where quantum chaos was seen to emerge (cf. Fig. 
[TJ. It also shows that (i) the predictions of the microcanonical 
ensemble become more accurate with increasing system size, 
and (ii) the increasing disagreement between the microcanon- 
ical prediction and the outcome of the relaxation dynamics, 
after crossing the superfluid to insulator transition, moves to 
larger values of V as the system size increases. 

To summarize, our studies indicate that thermalization does 
occur in gapped systems. If integrability is broken, ETH was 
shown to be valid, away from the edges of the spectrum, even 
if gaps are present in the spectrum and the ground state of 
the system is an insulator. We verified that ETH holds where 
quantum chaos develops and that thermalization closely fol- 
lows the validity of ETH (i.e., we found no instance where the 
rare event scenario put forward in Ref. lfl3ll emerges in these 
systems). Our analysis of different lattice sizes showed that: 
(i) The range of parameters over which ETH applies increases 
with increasing system size; in particular, ETH becomes valid 
deeper into the insulating side of the phase diagram, and (ii) 
away from integrability, the fluctuations of the eigenstate ex- 
pectation values of few-body observables decrease with in- 
creasing systems size. Further studies are needed to under- 
stand the dependence of the short-time dynamics on the initial 
state as well as the final Hamiltonian, and to determine the 
precise scaling, with system size, for the onset of ETH and 
quantum chaos. 
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I. HAMILTONIAN AND GAPLESS-GAPPED TRANSITION 



II. QUANTUM CHAOS INDICATORS 



We study the one-dimensional hard-core boson (HCB) 
model with nearest-neighbor (NN) hopping t, and nearest- 
and next-nearest-neighbor (NNN) interaction V and V', re- 
spectively. The Hamiltonian is written as 



(2) 



-V I - r 



i+l 2 



V n - i 



t+2 2 



Above, we take h — 1, L is the size of the chain, bi {b\ ) 
is the bosonic annihilation (creation) operator on site i and 
h\ = b\bi is the boson local density operator. HCBs do not 
occupy the same site, so b\ = b\ 2 = 0. 

Hamiltonian (f2]) conserves the total number of particles Nf, 
and is translational invariant; it is then composed of indepen- 
dent blocks each associated with a value of iVj and a total 
momentum k. Here we select Nf, = L/3 and study all values 
of k, from to ... \ L/2\. Lattices with up to 24 sites and 8 
HCBs are considered, corresponding to a total Hilbert space of 
dimension D — 735 471. The dimension D k of each fc-sector 
is given in Table J] We perform a full exact diagonalization of 
each sector independently. 



TABLE I: Dimension of subspaces 



L = 18 
D k 


k = 0,6 
1038 


k = 1,5,7 
1026 


fe = 2,4,8 
1035 


k = 3,9 
1028 


L = 21 
D k 


k = 0,7 
5538 


all other k's 
5537 






L = 24 
D k 


k = 0,8 
30667 


odd fc's 
30624 


k = 2,6,10 
30664 


k = 4 
30666 



In what follows, t = 1 sets the energy scale and we only 
consider repulsive interactions V, V > 0. We fix V = 6 and 
vary V' from to 9. The system is integrable when V = 0, 
while the addition of NNN interaction may lead to the onset 
of chaos. In addition, there is a critical value of the NNN 
interaction, V' c = 3, below which the ground state is a gapless 
superfluid and above which it becomes a gapped insulator (see 
Ref. [17] in main text). 



To analyze the transition from integrability to chaos, we 
study chaos indicators that depend only on the eigenvalues of 
the system, such as level spacing distribution and level number 
variance, and indicators that measure the level of derealiza- 
tion of the eigenvectors, such as the inverse participation ratio 
(IPR) and the information (Shannon) entropy (S). 



A. Level spacing distribution and level number variance 

Level spacing distribution and level number variance give 
information about short-range and long-range correlations, re- 
spectively. They are obtained from the unfolded spectrum of 
each symmetry sector separately. The procedure of unfold- 
ing consists of locally rescaling the energies, so that the mean 
level density of the new sequence of energies is 1 . 

For integrable systems, the distribution of spacings s of 
neighboring energy levels may cross and the distribution is 
Poissonian, 

P P (s) = cxp(-s), 

whereas in nonintegrable systems, level repulsion leads to the 
Wigner-Dyson distribution. The form of the latter depends 
on the symmetries of the system. Hamiltonian (O is time- 
reversal and rotationally invariant, it therefore gives the same 
distribution as random matrices with real and symmetric ele- 
ments, the so-called Gaussian Orthogonal Ensembles (GOEs): 



-Pgoe(s) 



■ exp 



7TS 



The level number variance is defined as 

X 2 (l)^((N(l,e) 2 )~(N(l,e))\ 

where N(l, e) is the number of states in the energy interval 
[e, e + 1} and (.) is the average over different initial values of 
the energy level e. For a Poisson distribution, 

£2,(0 = /, 
while for GOEs in the limit of large I, 



^GOe(0 — — 2 



ln(27r0 + 7+1 
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FIG. 5: (Color online) (a) Peak position of the level spacing distri- 
bution averaged over all sectors with k — 1, . . . , |_L/2J. The inset 
shows the quantity a (see text), (b) Level number variance averaged 
over the same fc-sectors, for L = 21, and compared to the Poisson 
and GOE results, (c) Shannon entropy in the mean-field (top) and 
momentum (bottom) basis vs energy per site; k — 2; black: L = 24, 
red: L = 21; dashed line: Sgoe ~ 9.5969 for D 2 = 30 664. 



where 7 is the Euler constant. 

Results for level spacing distribution and level number vari- 
ance for various values of V' are shown in Figs.|5ja) and|5jb). 
These figures are equivalent to Figs. 1(a) and 1(b) in the main 
text. Note that because parity is a symmetry found only in 
the sector with k = 0, this subspace is not included in the 
averages for both figures. 

In Fig. 1(a) of the main text we quantified the integrable- 
chaos transition in terms of the parameter /3 of the Brody dis- 
tribution. In Fig. |3 a) we consider the peak position of the 
distribution, which shifts from zero to w 0.7979 during the 
crossover from Pp(s) to Pwd{s), and a quantity a (inset) 
defined as 

a _ Ei\P^i)-PwD(8i)\ _ (3) 

Above the sum runs over the whole spectrum. In the chaotic 
limit a — > 0. 

FiguresQa) and|5jb) reinforce the main results of Figs. 1(a) 
and 1(b) in the main text, which are: (i) Two transitions are 
seen, from integrability to chaos as V' increases from V = 0, 
and a departure from chaoticity for large values of V . (ii) By 
comparing the values of V that induce chaos with the V for 
the superfluid-insulator transition in the inset of Fig. 1(a), a 
region of overlap between the gapped phase and the chaotic 
regime is identified, (iii) A direct dependence exists between 
the system size and the width of the interval of V values that 
lead to chaos; for larger systems the chaotic behavior goes 
deeper into the insulating phase. Larger system sizes bring 
also better agreement with the GOE results, as seen by com- 
paring Fig. 1(b) in the main text, which was obtained for 24 
sites, with Fig.|3b), which deals with L = 21. 



B. Derealization Measures 

Derealization measures, such as IPR and S, quantify the 
level of complexity of the eigenvectors (see references in 
main text). In general, they depend on the basis in which 
the computations are performed. For an eigenstate \ip a ) of 
Hamiltonian written in the basis vectors \(j>j) as \ip a ) = 

£j=i c q l^j)' an d S are respectively given by 
and 

D k 

S Q ^-£K| 2 lnK| 2 . 

i=i 

These quantities indicate how much spread each \ip a ) is in the 
selected basis. 

In the case of GOEs, the eigenvectors do not depend on the 
basis. They are simply random vectors, which give IPRgoe ~ 
Dfc/3 and Sgoe ~ ln(0.48£)fc). An essential difference of our 
system with respect to GOEs, besides the absence of random- 
ness on it, is that Hamiltonian <(2j has only two-body interac- 
tions. Consequently, its eigenstates may approach the GOE 
results only in the middle of the spectrum, close to the edges 
chaos does not fully develop. 

In Fig. HJc), we show S for two system sizes and in two 
bases. In the mean-field (mf) basis, \<frj)'s, correspond to 
the eigenstates of the integrable Hamiltonian (V = 0); this 
choice separates regular from chaotic behavior. In the fc-basis, 
\4>j)'s are the basis vectors of total momentum. As V' in- 
creases from zero and the system undergoes the first transi- 
tion, from integrability to chaos, the values of S m f become 
larger, indicating derealization of the eigenvectors in the mf- 
basis. The second transition, signaling the departure from 
chaoticity, is followed by the reduction of the values of Sfc 
and the consequent localization of the eigenvectors in the fc- 
basis. Even though the shrinking of the eigenvectors in fc- 
space is better visualized in terms of IPR (main text), it is 
important to observe again that, similarly to the distancing of 
the eigenvalues from Pwd{s) and S GOE (Z), this localization 
occurs when V is already beyond the critical point V c for the 
gapless-gapped transition. 

The structure of the eigenstates gives us information about 
what to expect in terms of thermalization. The eigenstate ther- 
malization hypothesis (ETH) states that thermalization should 
happen when the eigenstate expectation values (EEVs) do not 
fluctuate for states which are close in energy. In this case, 
the EEVs will equal the microcanonical average. The valid- 
ity of ETH is therefore certain to hold for GOEs, where all 
eigenstates are extended and have the same level of complex- 
ity. For our two-body interaction system, two aspects must 
be taken into account: chaoticity and the energy of the initial 
state. Away from the chaotic regime, that is when V' — > or 
V' 3> V c , the structure of eigenstates close in energy fluctu- 
ate significantly, as shown by the panels of Fig. [2c). But even 



7 



in the chaotic region, fluctuations are seen also at the edges 
of the spectrum. Thus, thermalization in our lattice should 
occur in the chaotic domain and for initial states with energy 
away from the spectrum borders. This conclusion is the same 
we arrived at in Ref. [16] (main text), where we studied sys- 
tems whose ground state was always gapless; therefore, as 
expected, the superfluid-insulator transition (a quantum phase 
transition) does not appear to affect the level of complexity of 
the bulk states. We note, however, that the results of the de- 
localization measures for system (f5]i show an asymmetry with 
respect to the middle of the spectrum that was not so evident in 
the systems of Ref. [16] (cf. Fig. 1(c), Fig.|5jc), and Figs. 1 1- 
16 from Ref. [16]). In the present lattice, larger fluctuations in 
the values of S and IPR are present for E > for all values of 
V studied. 



III. EIGENSTATE THERMALIZATION HYPOTHESIS 

The connection between the ETH and the chaos indicators 
is reinforced by comparing the outcomes for the derealiza- 
tion measures with those for the EEVs. 
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FIG. 6: (Color online) (a) Eigenstate expectation values (EEVs) of 
K (top) and n(k) (bottom) vs energy per site for the full spectrum 
(including all momentum sectors). Results are shown for four dif- 
ferent values of V' and L — 21. Panels (b) and (c) correspond to 
the average relative deviation of the EEVs of K (b) and n(k = 0) 
(c) with respect to the microcanonical result vs V' for T = 5 (see 
text) and three different lattice sizes. T — 5 (L — 24) corresponds 
to E = -6.24 for V' = 1, E = -5.55 for V = 3, E = -6.61 for 
V = 5, and E = -12.68 for V = 9. 

The top and bottom panels of Fig. |6ja) show, respectively, 
the EEVS of the kinetic energy 



K = J2~t [}th+i + H.c. 

i 

and the momentum distribution function 

m = jY, e ~ ikii ~ jrb l h 



for all eigenstates. Fig. |6ja) is analogous to Fig. 2(a), but 
refers to a smaller system. The results in both figures mir- 
ror the behavior of the derealization measures throughout 
the transitions. Large fluctuations of the EEVs are seen over 
the entire spectrum when the system is away from the chaotic 
limit, that is, when it approaches integrability and when it ap- 
proaches localization in fc-space. When chaos sets in, the fluc- 
tuations decrease mainly in the center of the spectrum, where 
the ETH is expected to be valid. We also should stress that, 
by comparing Fig. |6ja) with Fig. 2(a) in the chaotic regime, 
one can see that the fluctuations of the observables in all in- 
dividual eigenstates (away from the edges of the spectrum) 
decrease with increasing systems size. This is a clear indica- 
tion that the rare state scenario introduced in Ref. [13] in the 
main text does not take place in these systems. 

To quantify the deviation of the EEV for an observable O 
with respect to the microcanonical result (A"" c ), we define 



A mic O 



Xa Oaa 



Above, the sum runs over the microcanonical window, O aa 
are the EEVs of the operator O, and the microcanonical ex- 
pectation values Omic are computed as usual. 

Figures |6|b) and|6|c) show the relative deviation A mlc K 
and A mzc n(k = 0) averaged over all momentum sectors and 
for all eigenstates that lie within a window [E— AE, E+AE], 
where AE = 0.1. We select E according to the effective tem- 
perature T that we want to study [T = 5 here and T = 3 in the 
equivalent Figs. 2(b) and 2(c) in the main text]. The analysis 
in terms of a single temperature allows for a fair comparison 
of all systems sizes and values of V'. Temperature and energy 
are related by the expression 



(4) 



where Z = Tr |e #7 fc B T j is the partition function, and we 

set ks to unity. 

The results for A mic K and A mic n(k = 0) for any given 
systems size comply with the predictions of chaos measures. 
As the system size increases, the average deviations for both 
observables decrease and the width of the interval of values 
of V' for which the EEVs approach the thermal average in- 
creases. We should add that for some specific values of V', the 
relative deviations are seen to be larger in some larger system 
sizes, that is, the curves in Fig. [6] and Fig. 2 in the main text 
cross. As we show next, this is a finite size effect related to 
bands of eigenvalues that move, within the temperature scale 
considered in this work, as the system size is changed. 



IV. FLUCTUATION OF THE EEVS AS A FUNCTION OF 
TEMPERATURE 

The study of the fluctuation of the EEVs as a function of 
temperature gives further support to the discussions above. 

In Fig. |7J we present results for A mlc n(k = 0) vs T for 
nine different values of V' for systems with 24 and 21 lat- 
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FIG. 7: (Color online) Average relative deviation of eigenstate ex- 
pectation values with respect to the microcanonical prediction as a 
function of the effective temperature T of the eigenstates. Results 
for A mlc n(fc = 0) are given for the nine values of V' indicated and 
for systems with L = 24 (thick lines) and L — 21 (thin lines). 



tice sites and T < 10. As before, distinct features are as- 
sociated with different regimes. Far from chaoticity, large 
values of A mlc n(k = 0) appear for all temperatures consid- 
ered. In the chaotic domain, on the other hand, large values 
of A mlc n(k = 0) are restricted to low temperatures, while at 
large T, A mlc n(k = 0) saturates at small values. This corrob- 
orates our statements that the validity of ETH for these finite 
systems goes hand in hand with the onset of chaos and holds 
away from the edges of the spectrum. 

In general, we also observe that larger systems reduce the 
average relative deviations. However, especially away from 
chaoticity, this rule may not be obeyed for particular values of 
the effective temperature and particular values of V . Notice 
in Fig. [7] the peaks with large fluctuations that move within 
our temperature scale as the system size is increased. These 
are responsible for the bumps seen in Figs. 2(b) and 2(c) in 
the main text and for the behavior of the fluctuations of the 
system with L = 21 in Figs. [6fb) and|6jc). From the behav- 
ior seen here with increasing systems size, we expect them to 
disappear for very large system sizes. 



V. LONG-TIME DYNAMICS AFTER A QUANTUM 
QUENCH 

The endorsement of the predictions for thermalization 
drawn from chaos measures and the ETH is done in two steps. 
First we study the relaxation dynamics of the system to ver- 
ify that it brings the observables O close to the values of the 
diagonal ensemble, 



where C a is the overlap of the initial state 1^(0)) with the 
eigenstate a of the Hamiltonian, C a = (^^(O)). Then we 
compare the results for the diagonal ensemble with those for 
the microcanonical ensemble. When they coincide we say that 
thermalization has taken place. 

In the main text, the second step was presented in Fig. 4 
for two values of temperature, T = 3 and T = 5. The relative 
differences between the predictions of the microcanonical and 
diagonal ensembles for K and n(k) became indeed small for 
values of V' where the ETH was shown to be valid. The first 
step was illustrated in Fig. 3 for T = 3 and, for completeness, 
we now demonstrate that similar results hold for the dynamics 
of systems for which the effective temperature is T = 5. 
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FIG. 8: (Color online) Dynamics of the normalized difference be- 
tween the evolving expectation values of K (left panels) and n(k) 
(right panels) and the diagonal ensemble prediction. The plots are the 
result of the average over the evolution of nine initial states selected 
from the eigenstates of the Hamiltonian with Vl ni — 0, 1, . . . , 9 (ex- 
cluding the V used for the dynamics). The nine states for each 
V' were chosen such that the (conserved) energies during the time 
evolution are the same in all cases, and the effective temperature 
T = 5. Given the energy of the initial state in the final Hamil- 
tonian E = {ipini\H \ipini), the effective temperature is computed 
following Eq. 10. 



Figure [8] gives the normalized difference between the time 
evolving expectation value of K and n(k) and the diagonal 
ensemble prediction, 5 K (left panels) and Snk (right panels), 
respectively. We define SK(t) = \K(t) - K dlag \/\K diag \ 
and<5n fc (r) = (J2 k \ n ( k > T )- n diag(k)\)/(J2k n diag(k)). For 
each value of V considered in the analysis of the dynamics, 
nine initial states, all corresponding to T = 5, were selected 
from the eigenstates of the Hamiltonian with different values 
of V( ni (excluding V). We studied the evolution of the nine 
states after the quench and verified that their long-time dy- 
namics were very similar. In Fig. [8] we display the average 
over those nine different time evolutions. The plots show that 
after long times, observables relax to values similar to those 
predicted by the diagonal ensemble and those predictions be- 
come more accurate with increasing system size. Only for 
V = 9 we find large time fluctuations of K and a much 
longer time scale for relaxation, which is a consequence of 
the approach to localization in fc-space. However, even in the 
latter case, the time fluctuations are always seen to decrease 
with increasing system size. 



